home *** CD-ROM | disk | FTP | other *** search
/ CU Amiga Super CD-ROM 15 / CU Amiga Magazine's Super CD-ROM 15 (1997)(EMAP Images)(GB)[!][issue 1997-10].iso / CUCD / Graphics / Ghostscript / source / gxpcopy.c < prev    next >
C/C++ Source or Header  |  1997-06-10  |  24KB  |  807 lines
  1. /* Copyright (C) 1992, 1995, 1996, 1997 Aladdin Enterprises.  All rights reserved.
  2.   
  3.   This file is part of Aladdin Ghostscript.
  4.   
  5.   Aladdin Ghostscript is distributed with NO WARRANTY OF ANY KIND.  No author
  6.   or distributor accepts any responsibility for the consequences of using it,
  7.   or for whether it serves any particular purpose or works at all, unless he
  8.   or she says so in writing.  Refer to the Aladdin Ghostscript Free Public
  9.   License (the "License") for full details.
  10.   
  11.   Every copy of Aladdin Ghostscript must include a copy of the License,
  12.   normally in a plain ASCII text file named PUBLIC.  The License grants you
  13.   the right to copy, modify and redistribute Aladdin Ghostscript, but only
  14.   under certain conditions described in the License.  Among other things, the
  15.   License requires that the copyright notice and this notice be preserved on
  16.   all copies.
  17. */
  18.  
  19. /* gxpcopy.c */
  20. /* Path copying and flattening */
  21. #include "math_.h"
  22. #include "gx.h"
  23. #include "gserrors.h"
  24. #include "gconfigv.h"        /* for USE_FPU */
  25. #include "gxfixed.h"
  26. #include "gxfarith.h"
  27. #include "gzpath.h"
  28.  
  29. /* Forward declarations */
  30. private void adjust_point_to_tangent(P3(segment *, const segment *,
  31.                     const gs_fixed_point *));
  32. private int monotonize_internal(P2(gx_path *, const curve_segment *));
  33.  
  34. /* Copy a path, optionally flattening or monotonizing it. */
  35. /* If the copy fails, free the new path. */
  36. int
  37. gx_path_copy_reducing(const gx_path *ppath_old, gx_path *ppath,
  38.   fixed fixed_flatness, gx_path_copy_options options)
  39. {    gx_path old;
  40.     const segment *pseg;
  41.     int code;
  42.  
  43. #ifdef DEBUG
  44.     if ( gs_debug_c('P') )
  45.       gx_dump_path(ppath_old, "before reducing");
  46. #endif
  47.     old = *ppath_old;
  48.     if ( options & pco_init )
  49.       gx_path_init(ppath, ppath_old->memory);
  50.     pseg = (const segment *)(old.first_subpath);
  51.     while ( pseg )
  52.        {    switch ( pseg->type )
  53.           {
  54.           case s_start:
  55.             code = gx_path_add_point(ppath,
  56.                 pseg->pt.x, pseg->pt.y);
  57.             break;
  58.           case s_curve:
  59.            {    const curve_segment *pc = (const curve_segment *)pseg;
  60.             if ( fixed_flatness == max_fixed ) /* don't flatten */
  61.               { if ( options & pco_monotonize )
  62.                   code = monotonize_internal(ppath, pc);
  63.                 else
  64.                   code = gx_path_add_curve_notes(ppath,
  65.                 pc->p1.x, pc->p1.y, pc->p2.x, pc->p2.y,
  66.                 pc->pt.x, pc->pt.y, pseg->notes);
  67.               }
  68.             else
  69.               { fixed x0 = ppath->position.x;
  70.                 fixed y0 = ppath->position.y;
  71.                 int k = gx_curve_log2_samples(x0, y0, pc,
  72.                               fixed_flatness);
  73.                 segment_notes notes = pseg->notes;
  74.                 segment *start;
  75.                 curve_segment cseg;
  76.  
  77.                 if ( options & pco_accurate )
  78.                   { /* Add an extra line, which will become */
  79.                 /* the tangent segment. */
  80.                 code = gx_path_add_line_notes(ppath, x0, y0,
  81.                              notes);
  82.                 if ( code < 0 )
  83.                   break;
  84.                 start = ppath->current_subpath->last;
  85.                 notes |= sn_not_first;
  86.                   }
  87.                 cseg = *pc;
  88.                 code = gx_flatten_sample(ppath, k, &cseg, notes);
  89.                 if ( options & pco_accurate )
  90.                   { /*
  91.                  * Adjust the first and last segments so that
  92.                  * they line up with the tangents.
  93.                  */
  94.                 segment *end = ppath->current_subpath->last;
  95.  
  96.                 if ( code < 0 ||
  97.                      (code = gx_path_add_line_notes(ppath,
  98.                         ppath->position.x,
  99.                         ppath->position.y,
  100.                         pseg->notes | sn_not_first)) < 0
  101.                    )
  102.                   break;
  103.                 adjust_point_to_tangent(start, start->next,
  104.                             &pc->p1);
  105.                 adjust_point_to_tangent(end, end->prev,
  106.                             &pc->p2);
  107.                   }
  108.               }
  109.             break;
  110.            }
  111.           case s_line:
  112.             code = gx_path_add_line_notes(ppath,
  113.                 pseg->pt.x, pseg->pt.y, pseg->notes);
  114.             break;
  115.           case s_line_close:
  116.             code = gx_path_close_subpath(ppath);
  117.             break;
  118.           default:        /* can't happen */
  119.             code = gs_note_error(gs_error_unregistered);
  120.           }
  121.         if ( code < 0 )
  122.            {    gx_path_release(ppath);
  123.             if ( ppath == ppath_old )
  124.               *ppath = old;
  125.             return code;
  126.            }
  127.         pseg = pseg->next;
  128.     }
  129.     if ( path_last_is_moveto(&old) )
  130.       gx_path_add_point(ppath, old.position.x, old.position.y);
  131. #ifdef DEBUG
  132.     if ( gs_debug_c('P') )
  133.       gx_dump_path(ppath, "after reducing");
  134. #endif
  135.     return 0;
  136. }
  137.  
  138. /*
  139.  * Adjust one end of a line (the first or last line of a flattened curve)
  140.  * so it falls on the curve tangent.  The closest point on the line from
  141.  * (0,0) to (C,D) to a point (U,V) -- i.e., the point on the line at which
  142.  * a perpendicular line from the point intersects it -- is given by
  143.  *    T = (C*U + D*V) / (C^2 + D^2)
  144.  *    (X,Y) = (C*T,D*T)
  145.  * However, any smaller value of T will also work: the one we actually
  146.  * use is 0.25 * the value we just derived.  We must check that
  147.  * numerical instabilities don't lead to a negative value of T.
  148.  */
  149. private void
  150. adjust_point_to_tangent(segment *pseg, const segment *next,
  151.   const gs_fixed_point *p1)
  152. {    const fixed x0 = pseg->pt.x, y0 = pseg->pt.y;
  153.     const fixed fC = p1->x - x0, fD = p1->y - y0;
  154.  
  155.     /*
  156.      * By far the commonest case is that the end of the curve is
  157.      * horizontal or vertical.  Check for this specially, because
  158.      * we can handle it with far less work (and no floating point).
  159.      */
  160.     if ( fC == 0 ) {
  161.       /* Vertical tangent. */
  162.       const fixed DT = arith_rshift(next->pt.y - y0, 2);
  163.  
  164.       if ( fD == 0 )
  165.         return;        /* anomalous case */
  166.       if_debug1('2', "[2]adjusting vertical: DT = %g\n",
  167.             fixed2float(DT));
  168.       if ( DT > 0 )
  169.         pseg->pt.y = DT + y0;
  170.     } else if ( fD == 0 ) {
  171.       /* Horizontal tangent. */
  172.       const fixed CT = arith_rshift(next->pt.x - x0, 2);
  173.  
  174.       if_debug1('2', "[2]adjusting horizontal: CT = %g\n",
  175.             fixed2float(CT));
  176.       if ( CT > 0 )
  177.         pseg->pt.x = CT + x0;
  178.     } else {
  179.       /* General case. */
  180.       const double C = fC, D = fD;
  181.       const double T = (C * (next->pt.x - x0) + D * (next->pt.y - y0)) /
  182.         (C * C + D * D);
  183.  
  184.       if_debug3('2', "[2]adjusting: C = %g, D = %g, T = %g\n",
  185.             C, D, T);
  186.       if ( T > 0 ) {
  187.         pseg->pt.x = arith_rshift((fixed)(C * T), 2) + x0;
  188.         pseg->pt.y = arith_rshift((fixed)(D * T), 2) + y0;
  189.       }
  190.     }
  191. }
  192.  
  193. /* ---------------- Curve flattening ---------------- */
  194.  
  195. #define x1 pc->p1.x
  196. #define y1 pc->p1.y
  197. #define x2 pc->p2.x
  198. #define y2 pc->p2.y
  199. #define x3 pc->pt.x
  200. #define y3 pc->pt.y
  201.  
  202. #ifdef DEBUG
  203. private void
  204. dprint_curve(const char *str, fixed x0, fixed y0, const curve_segment *pc)
  205. {    dprintf9("%s p0=(%g,%g) p1=(%g,%g) p2=(%g,%g) p3=(%g,%g)\n",
  206.          str, fixed2float(x0), fixed2float(y0),
  207.          fixed2float(pc->p1.x), fixed2float(pc->p1.y),
  208.          fixed2float(pc->p2.x), fixed2float(pc->p2.y),
  209.          fixed2float(pc->pt.x), fixed2float(pc->pt.y));
  210. }
  211. #endif
  212.  
  213. /* Initialize a cursor for rasterizing a monotonic curve. */
  214. void
  215. gx_curve_cursor_init(curve_cursor *prc, fixed x0, fixed y0,
  216.   const curve_segment *pc, int k)
  217. {    fixed v01, v12;
  218.     int k2 = k + k, k3 = k2 + k;
  219.  
  220. #define bits_fit(v, n)\
  221.   (any_abs(v) <= max_fixed >> (n))
  222. /* The +2s are because of t3d and t2d, see below. */
  223. #define coeffs_fit(a, b, c)\
  224.   (k3 <= sizeof(fixed) * 8 - 3 &&\
  225.    bits_fit(a, k3 + 2) && bits_fit(b, k2 + 2) && bits_fit(c, k + 1))
  226.  
  227.     prc->k = k;
  228.     prc->p0.x = x0, prc->p0.y = y0;
  229.     prc->pc = pc;
  230.     /* Compute prc->a..c taking into account reversal of xy0/3 */
  231.     /* in curve_x_at_y. */
  232.     { fixed w0, w1, w2, w3;
  233.       if ( y0 < y3 )
  234.         w0 = x0, w1 = x1, w2 = x2, w3 = x3;
  235.       else
  236.         w0 = x3, w1 = x2, w2 = x1, w3 = x0;
  237.       curve_points_to_coefficients(w0, w1, w2, w3,
  238.                        prc->a, prc->b, prc->c, v01, v12);
  239.     }
  240.     prc->double_set = false;
  241.     prc->fixed_limit =
  242.       (coeffs_fit(prc->a, prc->b, prc->c) ? (1 << k) - 1 : -1);
  243.     /* Initialize the cache. */
  244.     prc->cache.ky0 = prc->cache.ky3 = y0;
  245.     prc->cache.xl = x0;
  246.     prc->cache.xd = 0;
  247. }
  248.  
  249. /*
  250.  * Determine the X value on a monotonic curve at a given Y value.  It turns
  251.  * out that we use so few points on the curve that it's actually faster to
  252.  * locate the desired point by recursive subdivision each time than to try
  253.  * to keep a cursor that we move incrementally.  What's even more surprising
  254.  * is that if floating point arithmetic is reasonably fast, it's faster to
  255.  * compute the X value at the desired point explicitly than to do the
  256.  * recursive subdivision on X as well as Y.
  257.  */
  258. #define SUBDIVIDE_X USE_FPU_FIXED
  259. fixed
  260. gx_curve_x_at_y(curve_cursor *prc, fixed y)
  261. {    fixed xl, xd;
  262.     fixed yd, yrel;
  263.  
  264.     /* Check the cache before doing anything else. */
  265.     if ( y >= prc->cache.ky0 && y <= prc->cache.ky3 )
  266.       { yd = prc->cache.ky3 - prc->cache.ky0;
  267.         yrel = y - prc->cache.ky0;
  268.         xl = prc->cache.xl;
  269.         xd = prc->cache.xd;
  270.         goto done;
  271.       }
  272.  
  273.     {
  274. #define x0 prc->p0.x
  275. #define y0 prc->p0.y
  276.     const curve_segment *pc = prc->pc;
  277.     /* Reduce case testing by ensuring y3 >= y0. */
  278.     fixed cy0 = y0, cy1, cy2, cy3 = y3;
  279.     fixed cx0;
  280. #if SUBDIVIDE_X
  281.     fixed cx1, cx2, cx3;
  282. #else
  283.     int t = 0;
  284. #endif
  285.     int k, i;
  286.  
  287.     if ( cy0 > cy3 )
  288.       cx0 = x3,
  289. #if SUBDIVIDE_X
  290.         cx1 = x2, cx2 = x1, cx3 = x0,
  291. #endif
  292.         cy0 = y3, cy1 = y2, cy2 = y1, cy3 = y0;
  293.     else
  294.       cx0 = x0,
  295. #if SUBDIVIDE_X
  296.         cx1 = x1, cx2 = x2, cx3 = x3,
  297. #endif
  298.         cy1 = y1, cy2 = y2;
  299. #define midpoint_fast(a,b)\
  300.   arith_rshift_1((a) + (b) + 1)
  301.     for ( i = k = prc->k; i > 0; --i )
  302.       {
  303.         fixed ym = midpoint_fast(cy1, cy2);
  304.         fixed yn = ym + arith_rshift(cy0 - cy1 - cy2 + cy3 + 4, 3);
  305. #if SUBDIVIDE_X
  306.         fixed xm = midpoint_fast(cx1, cx2);
  307.         fixed xn = xm + arith_rshift(cx0 - cx1 - cx2 + cx3 + 4, 3);
  308. #else
  309.         t <<= 1;
  310. #endif
  311.  
  312.         if ( y < yn )
  313. #if SUBDIVIDE_X
  314.           cx1 = midpoint_fast(cx0, cx1),
  315.           cx2 = midpoint_fast(cx1, xm),
  316.           cx3 = xn,
  317. #endif
  318.           cy1 = midpoint_fast(cy0, cy1),
  319.           cy2 = midpoint_fast(cy1, ym),
  320.           cy3 = yn;
  321.         else
  322. #if SUBDIVIDE_X
  323.           cx2 = midpoint_fast(cx2, cx3),
  324.           cx1 = midpoint_fast(xm, cx2),
  325.           cx0 = xn,
  326. #else
  327.           t++,
  328. #endif
  329.           cy2 = midpoint_fast(cy2, cy3),
  330.           cy1 = midpoint_fast(ym, cy2),
  331.           cy0 = yn;
  332.       }
  333. #if SUBDIVIDE_X
  334.     xl = cx0;
  335.     xd = cx3 - cx0;
  336. #else
  337.     { fixed a = prc->a, b = prc->b, c = prc->c;
  338.  
  339. /* We must use (1 << k) >> 1 instead of 1 << (k - 1) in case k == 0. */
  340. #define compute_fixed(a, b, c)\
  341.   arith_rshift(arith_rshift(arith_rshift(a * t3, k) + b * t2, k)\
  342.            + c * t + ((1 << k) >> 1), k)
  343. #define compute_diff_fixed(a, b, c)\
  344.   arith_rshift(arith_rshift(arith_rshift(a * t3d, k) + b * t2d, k)\
  345.            + c, k)
  346.  
  347.                 /* use multiply if possible */
  348. #define np2(n) (1.0 / (1L << (n)))
  349.           static const double k_denom[11] =
  350.         { np2(0), np2(1), np2(2), np2(3), np2(4),
  351.           np2(5), np2(6), np2(7), np2(8), np2(9), np2(10)
  352.         };
  353.           static const double k2_denom[11] =
  354.         { np2(0), np2(2), np2(4), np2(6), np2(8),
  355.           np2(10), np2(12), np2(14), np2(16), np2(18), np2(20)
  356.         };
  357.           static const double k3_denom[11] =
  358.         { np2(0), np2(3), np2(6), np2(9), np2(12),
  359.           np2(15), np2(18), np2(21), np2(24), np2(27), np2(30)
  360.         };
  361.           double den1, den2;
  362. #undef np2
  363.  
  364. #define setup_floating(da, db, dc, a, b, c)\
  365.   (k >= countof(k_denom) ?\
  366.    (den1 = ldexp(1.0, -k),\
  367.     den2 = den1 * den1,\
  368.     da = (den2 * den1) * a,\
  369.     db = den2 * b,\
  370.     dc = den1 * c) :\
  371.    (da = k3_denom[k] * a,\
  372.     db = k2_denom[k] * b,\
  373.     dc = k_denom[k] * c))
  374. #define compute_floating(da, db, dc)\
  375.   ((fixed)(da * t3 + db * t2 + dc * t + 0.5))
  376. #define compute_diff_floating(da, db, dc)\
  377.   ((fixed)(da * t3d + db * t2d + dc))
  378.  
  379.       if ( t <= prc->fixed_limit )
  380.         { /* We can do everything in integer/fixed point. */
  381.           int t2 = t * t, t3 = t2 * t;
  382.           int t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  383.           xl = compute_fixed(a, b, c) + cx0;
  384.           xd = compute_diff_fixed(a, b, c);
  385. #ifdef DEBUG
  386.           { double fa, fb, fc;
  387.         fixed xlf, xdf;
  388.         setup_floating(fa, fb, fc, a, b, c);
  389.         xlf = compute_floating(fa, fb, fc) + cx0;
  390.         xdf = compute_diff_floating(fa, fb, fc);
  391.         if ( any_abs(xlf - xl) > fixed_epsilon ||
  392.              any_abs(xdf - xd) > fixed_epsilon
  393.            )
  394.           dprintf9("Curve points differ: k=%d t=%d a,b,c=%g,%g,%g\n   xl,xd fixed=%g,%g floating=%g,%g\n",
  395.                k, t,
  396.                fixed2float(a), fixed2float(b), fixed2float(c),
  397.                fixed2float(xl), fixed2float(xd),
  398.                fixed2float(xlf), fixed2float(xdf));
  399. /*xl = xlf, xd = xdf;*/
  400.           }
  401. #endif
  402.         }
  403.       else
  404.         { /*
  405.            * Either t3 (and maybe t2) won't fit in an int, or more
  406.            * likely the result of the multiplies won't fit.
  407.            */
  408. #define fa prc->da
  409. #define fb prc->db
  410. #define fc prc->dc
  411.           if ( !prc->double_set )
  412.         { setup_floating(fa, fb, fc, a, b, c);
  413.           prc->double_set = true;
  414.         }
  415.           if ( t < 1L << ((sizeof(long) * 8 - 1) / 3) )
  416.         { /*
  417.            * t3 (and maybe t2) might not fit in an int, but they
  418.            * will fit in a long.  If we have slow floating point,
  419.            * do the computation in double-precision fixed point,
  420.            * otherwise do it in fixed point.
  421.            */
  422.           long t2 = (long)t * t, t3 = t2 * t;
  423.           long t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  424.  
  425.           xl = compute_floating(fa, fb, fc) + cx0;
  426.           xd = compute_diff_floating(fa, fb, fc);
  427.         }
  428.           else
  429.         { /*
  430.            * t3 (and maybe t2) don't even fit in a long.
  431.            * Do the entire computation in floating point.
  432.            */
  433.           double t2 = (double)t * t, t3 = t2 * t;
  434.           double t3d = (t2 + t) * 3 + 1, t2d = t + t + 1;
  435.  
  436.           xl = compute_floating(fa, fb, fc) + cx0;
  437.           xd = compute_diff_floating(fa, fb, fc);
  438.         }
  439. #undef fa
  440. #undef fb
  441. #undef fc
  442.         }
  443.     }
  444. #endif                /* (!)SUBDIVIDE_X */
  445.  
  446.     /* Update the cache. */
  447.     prc->cache.ky0 = cy0;
  448.     prc->cache.ky3 = cy3;
  449.     prc->cache.xl = xl;
  450.     prc->cache.xd = xd;
  451.     yd = cy3 - cy0;
  452.     yrel = y - cy0;
  453. #undef x0
  454. #undef y0
  455.     }
  456. done:
  457.     /*
  458.      * Now interpolate linearly between current and next.
  459.      * We know that 0 <= yrel < yd.
  460.      * It's unlikely but possible that cy0 = y = cy3:
  461.      * handle this case specially.
  462.      */
  463.     if ( yrel == 0 )
  464.       return xl;
  465.     /*
  466.      * Compute in fixed point if possible.
  467.      */
  468. #define half_fixed_bits ((fixed)1 << (sizeof(fixed) * 4))
  469.     if ( yrel < half_fixed_bits )
  470.       { if ( xd >= 0 )
  471.           { if ( xd < half_fixed_bits )
  472.           return (ufixed)xd * (ufixed)yrel / (ufixed)yd + xl;
  473.           }
  474.         else
  475.           { if ( xd > -half_fixed_bits )
  476.           return -(fixed)((ufixed)(-xd) * (ufixed)yrel / (ufixed)yd) + xl;
  477.           }
  478.       }
  479. #undef half_fixed_bits
  480.     return fixed_mult_quo(xd, yrel, yd) + xl;
  481. }
  482.  
  483. #undef x1
  484. #undef y1
  485. #undef x2
  486. #undef y2
  487. #undef x3
  488. #undef y3
  489.  
  490. /* ---------------- Monotonic curves ---------------- */
  491.  
  492. /* Test whether a path is free of non-monotonic curves. */
  493. bool
  494. gx_path_is_monotonic(const gx_path *ppath)
  495. {    const segment *pseg = (const segment *)(ppath->first_subpath);
  496.     gs_fixed_point pt0;
  497.  
  498.     while ( pseg )
  499.        {    switch ( pseg->type )
  500.           {
  501.           case s_start:
  502.             {    const subpath *psub = (const subpath *)pseg;
  503.             /* Skip subpaths without curves. */
  504.             if ( !psub->curve_count )
  505.               pseg = psub->last;
  506.             }
  507.             break;
  508.           case s_curve:
  509.             {    const curve_segment *pc = (const curve_segment *)pseg;
  510.             double t[2];
  511.             int nz = gx_curve_monotonic_points(pt0.y,
  512.                     pc->p1.y, pc->p2.y, pc->pt.y, t);
  513.             if ( nz != 0 )
  514.               return false;
  515.             nz = gx_curve_monotonic_points(pt0.x,
  516.                     pc->p1.x, pc->p2.x, pc->pt.x, t);
  517.             if ( nz != 0 )
  518.               return false;
  519.             }
  520.             break;
  521.           default:
  522.             ;
  523.           }
  524.         pt0 = pseg->pt;
  525.         pseg = pseg->next;
  526.     }
  527.     return true;
  528. }
  529.  
  530. /* Monotonize a curve, by splitting it if necessary. */
  531. /* In the worst case, this could split the curve into 9 pieces. */
  532. private int
  533. monotonize_internal(gx_path *ppath, const curve_segment *pc)
  534. {    fixed x0 = ppath->position.x, y0 = ppath->position.y;
  535.     segment_notes notes = pc->notes;
  536.     double t[2];
  537. #define max_segs 9
  538.     curve_segment cs[max_segs];
  539.     const curve_segment *pcs;
  540.     curve_segment *pcd;
  541.     int i, j, nseg;
  542.     int nz;
  543.  
  544.     /* Monotonize in Y. */
  545.     nz = gx_curve_monotonic_points(y0, pc->p1.y, pc->p2.y, pc->pt.y, t);
  546.     nseg = max_segs - 1 - nz;
  547.     pcd = cs + nseg;
  548.     if ( nz == 0 )
  549.       *pcd = *pc;
  550.     else
  551.       { gx_curve_split(x0, y0, pc, t[0], pcd, pcd + 1);
  552.         if ( nz == 2 )
  553.           gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
  554.                  (t[1] - t[0]) / (1 - t[0]),
  555.                  pcd + 1, pcd + 2);
  556.       }
  557.  
  558.     /* Monotonize in X. */
  559.     for ( pcs = pcd, pcd = cs, j = nseg; j < max_segs; ++pcs, ++j )
  560.       { nz = gx_curve_monotonic_points(x0, pcs->p1.x, pcs->p2.x,
  561.                        pcs->pt.x, t);
  562.         
  563.         if ( nz == 0 )
  564.           *pcd = *pcs;
  565.         else
  566.           { gx_curve_split(x0, y0, pcs, t[0], pcd, pcd + 1);
  567.         if ( nz == 2 )
  568.           gx_curve_split(pcd->pt.x, pcd->pt.y, pcd + 1,
  569.                  (t[1] - t[0]) / (1 - t[0]),
  570.                  pcd + 1, pcd + 2);
  571.           }
  572.         pcd += nz + 1;
  573.         x0 = pcd[-1].pt.x;
  574.         y0 = pcd[-1].pt.y;
  575.       }
  576.     nseg = pcd - cs;
  577.  
  578.     /* Add the segment(s) to the output. */
  579. #ifdef DEBUG
  580.     if ( gs_debug_c('2') )
  581.       { int pi;
  582.         gs_fixed_point pp0;
  583.  
  584.         pp0 = ppath->position;
  585.         if ( nseg == 1 )
  586.           dprint_curve("[2]No split", pp0.x, pp0.y, pc);
  587.         else
  588.           { dprintf1("[2]Split into %d segments:\n", nseg);
  589.         dprint_curve("[2]Original", pp0.x, pp0.y, pc);
  590.         for ( pi = 0; pi < nseg; ++pi )
  591.           { dprint_curve("[2] =>", pp0.x, pp0.y, cs + pi);
  592.             pp0 = cs[pi].pt;
  593.           }
  594.           }
  595.       }
  596. #endif
  597.     for ( pcs = cs, i = 0; i < nseg; ++pcs, ++i )
  598.       { int code = gx_path_add_curve_notes(ppath, pcs->p1.x, pcs->p1.y,
  599.                            pcs->p2.x, pcs->p2.y,
  600.                            pcs->pt.x, pcs->pt.y,
  601.                            notes |
  602.                            (i > 0 ? sn_not_first :
  603.                         sn_none));
  604.         if ( code < 0 )
  605.           return code;
  606.       }
  607.  
  608.     return 0;
  609. }
  610.  
  611. /*
  612.  * Split a curve if necessary into pieces that are monotonic in X or Y as a
  613.  * function of the curve parameter t.  This allows us to rasterize curves
  614.  * directly without pre-flattening.  This takes a fair amount of analysis....
  615.  * Store the values of t of the split points in pst[0] and pst[1].  Return
  616.  * the number of split points (0, 1, or 2).
  617.  */
  618. int
  619. gx_curve_monotonic_points(fixed v0, fixed v1, fixed v2, fixed v3,
  620.   double pst[2])
  621. {
  622.     /*
  623.        Let
  624.         v(t) = a*t^3 + b*t^2 + c*t + d, 0 <= t <= 1.
  625.        Then
  626.         dv(t) = 3*a*t^2 + 2*b*t + c.
  627.        v(t) has a local minimum or maximum (or inflection point)
  628.        precisely where dv(t) = 0.  Now the roots of dv(t) = 0 (i.e.,
  629.        the zeros of dv(t)) are at
  630.           t =    ( -2*b +/- sqrt(4*b^2 - 12*a*c) ) / 6*a
  631.             =    ( -b +/- sqrt(b^2 - 3*a*c) ) / 3*a
  632.        (Note that real roots exist iff b^2 >= 3*a*c.)
  633.        We want to know if these lie in the range (0..1).
  634.        (The endpoints don't count.)  Call such a root a "valid zero."
  635.        Since computing the roots is expensive, we would like to have
  636.        some cheap tests to filter out cases where they don't exist
  637.        (i.e., where the curve is already monotonic).
  638.      */
  639.     fixed v01, v12, a, b, c, b2, a3;
  640.     fixed dv_end, b2abs, a3abs;
  641.  
  642.     curve_points_to_coefficients(v0, v1, v2, v3, a, b, c, v01, v12);
  643.     b2 = b << 1;
  644.     a3 = (a << 1) + a;
  645.     /*
  646.        If a = 0, the only possible zero is t = -c / 2*b.
  647.        This zero is valid iff sign(c) != sign(b) and 0 < |c| < 2*|b|.
  648.      */
  649.     if ( a == 0 )
  650.       { if ( (b ^ c) < 0 && any_abs(c) < any_abs(b2) && c != 0 )
  651.           { *pst = (double)(-c) / b2;
  652.         return 1;
  653.           }
  654.         else
  655.           return 0;
  656.       }
  657.     /*
  658.        Iff a curve is horizontal at t = 0, c = 0.  In this case,
  659.        there can be at most one other zero, at -2*b / 3*a.
  660.        This zero is valid iff sign(a) != sign(b) and 0 < 2*|b| < 3*|a|.
  661.      */
  662.     if ( c == 0 )
  663.       { if ( (a ^ b) < 0 && any_abs(b2) < any_abs(a3) && b != 0 )
  664.           { *pst = (double)(-b2) / a3;
  665.         return 1;
  666.           }
  667.         else
  668.           return 0;
  669.       }
  670.     /*
  671.        Similarly, iff a curve is horizontal at t = 1, 3*a + 2*b + c = 0.
  672.        In this case, there can be at most one other zero,
  673.        at -1 - 2*b / 3*a, iff sign(a) != sign(b) and 1 < -2*b / 3*a < 2,
  674.        i.e., 3*|a| < 2*|b| < 6*|a|.
  675.      */
  676.     else if ( (dv_end = a3 + b2 + c) == 0 )
  677.       { if ( (a ^ b) < 0 &&
  678.          (b2abs = any_abs(b2)) > (a3abs = any_abs(a3)) &&
  679.          b2abs < a3abs << 1
  680.            )
  681.           { *pst = (double)(-b2 - a3) / a3;
  682.         return 1;
  683.           }
  684.         else
  685.           return 0;
  686.       }
  687.     /*
  688.        If sign(dv_end) != sign(c), at least one valid zero exists,
  689.        since dv(0) and dv(1) have opposite signs and hence
  690.        dv(t) must be zero somewhere in the interval [0..1].
  691.      */
  692.     else if ( (dv_end ^ c) < 0 )
  693.       ;
  694.     /*
  695.        If sign(a) = sign(b), no valid zero exists,
  696.        since dv is monotonic on [0..1] and has the same sign
  697.        at both endpoints.
  698.      */
  699.     else if ( (a ^ b) >= 0 )
  700.       return 0;
  701.     /*
  702.        Otherwise, dv(t) may be non-monotonic on [0..1]; it has valid zeros
  703.        iff its sign anywhere in this interval is different from its sign
  704.        at the endpoints, which occurs iff it has an extremum in this
  705.        interval and the extremum is of the opposite sign from c.
  706.        To find this out, we look for the local extremum of dv(t)
  707.        by observing
  708.         d2v(t) = 6*a*t + 2*b
  709.        which has a zero only at
  710.         t1 = -b / 3*a
  711.        Now if t1 <= 0 or t1 >= 1, no valid zero exists.
  712.        Note that we just determined that sign(a) != sign(b), so we know t1 > 0.
  713.      */
  714.     else if ( any_abs(b) >= any_abs(a3) )
  715.       return 0;
  716.     /*
  717.        Otherwise, we just go ahead with the computation of the roots,
  718.        and test them for being in the correct range.  Note that a valid
  719.        zero is an inflection point of v(t) iff d2v(t) = 0; we don't
  720.        bother to check for this case, since it's rare.
  721.      */
  722.     { double nbf = (double)(-b);
  723.       double a3f = (double)a3;
  724.       double radicand = nbf * nbf - a3f * c;
  725.  
  726.       if ( radicand < 0 )
  727.         { if_debug1('2', "[2]negative radicand = %g\n", radicand);
  728.           return 0;
  729.         }
  730.       { double root = sqrt(radicand);
  731.         int nzeros = 0;
  732.         double z = (nbf - root) / a3f;
  733.  
  734.         /*
  735.          * We need to return the zeros in the correct order.
  736.          * We know that root is non-negative, but a3f may be either
  737.          * positive or negative, so we need to check the ordering
  738.          * explicitly.
  739.          */
  740.         if_debug2('2', "[2]zeros at %g, %g\n", z, (nbf + root) / a3f);
  741.         if ( z > 0 && z < 1 )
  742.           *pst = z, nzeros = 1;
  743.         if ( root != 0 )
  744.           { z = (nbf + root) / a3f;
  745.         if ( z > 0 && z < 1 )
  746.           { if ( nzeros && a3f < 0 ) /* order is reversed */
  747.               pst[1] = *pst, *pst = z;
  748.             else
  749.               pst[nzeros] = z;
  750.             nzeros++;
  751.           }
  752.           }
  753.         return nzeros;
  754.       }
  755.     }
  756. }
  757.  
  758. /*
  759.  * Split a curve at an arbitrary point t.  The above midpoint split is a
  760.  * special case of this with t = 0.5.
  761.  */
  762. void
  763. gx_curve_split(fixed x0, fixed y0, const curve_segment *pc, double t,
  764.   curve_segment *pc1, curve_segment *pc2)
  765. {    /*
  766.      * If the original function was v(t), we want to compute the points
  767.      * for the functions v1(T) = v(t * T) and v2(T) = v(t + (1 - t) * T).
  768.      * Straightforwardly,
  769.      *    v1(T) = a*t^3*T^3 + b*t^2*T^2 + c*t*T + d
  770.      * i.e.
  771.      *    a1 = a*t^3, b1 = b*t^2, c1 = c*t, d1 = d.
  772.      * Similarly,
  773.      *    v2(T) = a*[t + (1-t)*T]^3 + b*[t + (1-t)*T]^2 +
  774.      *        c*[t + (1-t)*T] + d
  775.      *          = a*[(1-t)^3*T^3 + 3*t*(1-t)^2*T^2 + 3*t^2*(1-t)*T +
  776.      *           t^3] + b*[(1-t)^2*T^2 + 2*t*(1-t)*T + t^2] +
  777.      *           c*[(1-t)*T + t] + d
  778.      *          = a*(1-t)^3*T^3 + [a*3*t + b]*(1-t)^2*T^2 +
  779.      *           [a*3*t^2 + b*2*t + c]*(1-t)*T +
  780.      *           a*t^3 + b*t^2 + c*t + d
  781.      * We do this in the simplest way, namely, we convert the points to
  782.      * coefficients, do the arithmetic, and convert back.  It would
  783.      * obviously be faster to do the arithmetic directly on the points,
  784.      * as the midpoint code does; this is just an implementation issue
  785.      * that we can revisit if necessary.
  786.      */
  787.     double t2 = t * t, t3 = t2 * t;
  788.     double omt = 1 - t, omt2 = omt * omt, omt3 = omt2 * omt;
  789.     fixed v01, v12, a, b, c, na, nb, nc;
  790.  
  791.     if_debug1('2', "[2]splitting at t = %g\n", t);
  792. #define compute_seg(v0, v)\
  793.     curve_points_to_coefficients(v0, pc->p1.v, pc->p2.v, pc->pt.v,\
  794.                      a, b, c, v01, v12);\
  795.     na = (fixed)(a * t3), nb = (fixed)(b * t2), nc = (fixed)(c * t);\
  796.     curve_coefficients_to_points(na, nb, nc, v0,\
  797.                      pc1->p1.v, pc1->p2.v, pc1->pt.v);\
  798.     na = (fixed)(a * omt3);\
  799.     nb = (fixed)((a * t * 3 + b) * omt2);\
  800.     nc = (fixed)((a * t2 * 3 + b * 2 * t + c) * omt);\
  801.     curve_coefficients_to_points(na, nb, nc, pc1->pt.v,\
  802.                      pc2->p1.v, pc2->p2.v, pc2->pt.v)
  803.     compute_seg(x0, x);
  804.     compute_seg(y0, y);
  805. #undef compute_seg
  806. }
  807.